COMMENT ⊗   VALID 00013 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00003 00002	\input acphdr % Section 3.2
C00004 00003	%folio 9 galley 4 (C) Addison-Wesley 1978	*
C00017 00004	%folio 13 galley 5 (C) Addison-Wesley 1978	*
C00024 00005	%folio 15 galley 6 (C) Addison-Wesley 1978	*
C00036 00006	%folio 20 galley 7 (C) Addison-Wesley 1978	*
C00049 00007	%folio 23 galley 8 (C) Addison-Wesley 1978	*
C00065 00008	%folio 26 galley 9 (C) Addison-Wesley 1978	*
C00079 00009	%folio 32 galley 10 (C) Addison-Wesley 1978	*
C00089 00010	%folio 37 galley 11 (C) Addison-Wesley 1978	*
C00101 00011	%folio 42 galley 12 (C) Addison-Wesley 1978	*
C00120 00012	%folio 46 galley 13 (C) Addison-Wesley 1978	*
C00139 00013	\vfill\eject
C00140 ENDMK
C⊗;
\input acphdr % Section 3.2
\runninglefthead{RANDOM NUMBERS} % chapter title
\titlepage\setcount00
\null
\vfill
\tenpoint
\ctrline{SECTION 3.2 of THE ART OF COMPUTER PROGRAMMING}
\ctrline{$\copyright$ 1978 Addison--Wesley Publishing Company, Inc.}
\vfill
\runningrighthead{GENERATING UNIFORM RANDOM NUMBERS}
\section{3.2}
\eject
\setcount0 9
%folio 9 galley 4 (C) Addison-Wesley 1978	*
\sectionbegin{3.2. GENERATING UNIFORM RANDOM NUMBERS}
I{\:cN THIS SECTION} we shall consider methods for
generating a sequence of random fractions, i.e., random {\sl
real numbers $U↓n$, uniformly distributed between zero and one.}
Since a computer can represent a real number with only finite
accuracy, we shall actually be generating integers $X↓n$ between
zero and some number $m$; the fraction
$$U↓n = X↓n/m$$ will then lie between zero
and one. Usually $m$ is the word size of the computer, so $X↓n$
may be regarded (conservatively) as the integer contents of
a computer word with the radix point assumed at the extreme
right, and $U↓n$ may be regarded (liberally) as the contents
of the same word with the radix point assumed at the extreme
left.

\runningrighthead{THE LINEAR CONGRUENTIAL METHOD} \section{3.2.1}
\sectionskip
\sectionbegin{3.2.1. The Linear Congruential Method}
By far the most popular random-number generators
in use today are special cases of the following scheme, introduced
by D. H. Lehmer in 1949.\xskip [See {\sl Proc.\ 2nd Symp.\ on Large-Scale
Digital Calculating Machinery} (Cambridge: Harvard University
Press, 1951), 141--146.]\xskip We choose four ``magic numbers'':
$$\vcenter{\halign{\rt{$#$}⊗\lft{\quad#}⊗\quad\rt{$#$}⊗\lft{$\;#$}\cr
m,⊗the modulus;⊗m⊗>0.\cr
a,⊗the multiplier;⊗0⊗≤a<m.\cr
c,⊗the increment;⊗0⊗≤c<m.\cr
X↓0,⊗the starting value;⊗0⊗≤X↓0<m.\cr
}} \eqno(1)$$
The desired sequence of random numbers $\langle X↓n \rangle$
is then obtained by setting
$$X↓{n+1} = (aX↓n + c)\mod m,\qquad n ≥ 0.\eqno (2)$$
This is called a {\sl linear congruential sequence.}
Taking the remainder mod $m$ is somewhat like determining where
a ball will land in a spinning roulette wheel.

For example, the sequence obtained when $m = 10$, $X↓0 = a = c
= 7$, is
$$7,\ 6,\ 9,\ 0,\ 7,\ 6,\ 9,\ 0,\ \ldotss\,. \eqno (3)$$ As this
example shows, the sequence is not always ``random'' for all
choices of $m$, $a$, $c$, and $X↓0$; the principles of choosing
the magic numbers appropriately will be investigated carefully in
later parts of this chapter.

Example (3) illustrates the fact that the congruential sequences
always ``get into a loop''; i.e., there is ultimately a cycle
of numbers that is repeated endlessly. This property is common
to all sequences having the general form $X↓{n+1} = f(X↓n)$;
see exercise 3.1--6. The repeating cycle is called the {\sl
period}\/; sequence (3) has a period of length 4. A useful
sequence will of course have a relatively long period.

The special case $c = 0$ deserves explicit mention, since the number generation
process is a little faster when $c = 0$ than it is when $c ≠
0$. We shall see later that the restriction $c = 0$ cuts
down the length of the period of the sequence, but it is still
possible to make the period reasonably long. Lehmer's original
generation method had $c = 0$, although he mentioned $c ≠ 0$
as a possibility; the idea of taking $c ≠ 0$
to obtain longer periods is due to Thomson [{\sl Comp.\ J. \bf 1}
(1958), 83, 86] and, independently, to Rotenberg [{\sl JACM \bf 7}
(1960), 75--77]. The terms {\sl multiplicative congruential
method} and {\sl mixed congruential method} are used by many
authors to denote linear congruential sequences with $c = 0$ and
$c ≠ 0$, respectively.

The letters $m$, $a$, $c$, and $X↓0$ will be used throughout this
chapter in the sense described above. Furthermore, we will find
it useful to define
$$b = a - 1,\eqno (4)$$ in order to simplify many
of our formulas.

We can immediately reject the case $a = 1$, for this would mean
that $X↓n = (X↓0 + nc)\mod m$, and the sequence would certainly
not behave as a random sequence. The case $a = 0$ is even worse.
Hence for practical purposes we may assume that
$$a ≥ 2,\qquad b ≥ 1.\eqno (5)$$
Now we can prove a generalization of Eq.\ (2),
$$X↓{n+k} = \left( a↑kX↓n + (a↑k - 1)c/b\right)\mod m,
\qquad k ≥ 0,\quad n ≥ 0,\eqno (6)$$
which expresses the $(n + k)$th term directly in terms
of the $n$th term.\xskip (The special case $n = 0$ in this equation
is worthy of note.)\xskip
It follows that the sub\-sequence consisting of every $k$th
term of our sequence is another linear congruential sequence,
having the multiplier $a↑k \mod m$ and the increment 
$\biglp(a↑k - 1)c/b\bigrp\mod m$.

An important corollary of (6) is that the general sequence defined
by $m$, $a$, $c$, and $X↓0$ can be expressed very simply in terms
of the special case where $c = 1$ and $X↓0 = 0$. Let
$$Y↓0 = 0,\qquad Y↓{n+1} = (aY↓n + 1) \mod m.\eqno (7)$$
According to Eq.\ (6) we will have $Y↓k=(a↑k-1)/b$, hence
the general sequence (2) satisfies
$$X↓n = (AY↓n + X↓0) \mod m,\qquad \hbox{where }
A = (X↓0b + c)\mod m.\eqno (8)$$ 

\exbegin{EXERCISES}
\exno 1. [10] Example
(3) shows a situation in which $X↓4 = X↓0$, so the sequence
begins again from the beginning. Give an example of a linear
congruential sequence with $m = 10$ for which $X↓0$ never appears
again in the sequence.

\trexno 2. [M20] Show that if $a$ and $m$ are relatively prime,
the number $X↓0$ will always appear in the period.

\exno 3. [M10] If $a$ and $m$ are not relatively prime, explain
why the sequence will be somewhat handicapped and probably not
very random; hence we will generally want the multiplier $a$ to be relatively
prime to the modulus $m$.

\exno 4. [11] Prove Eq.\ (6).

\exno 5. [M20] Equation
(6) holds for $k ≥ 0$. If possible, give a formula that expresses
$X↓{n+k}$ in terms of $X↓n$ for {\sl negative} values of $k$.

\runningrighthead{CHOICE OF MODULUS} \section{3.2.1.1}
\dimsectionbegin{3.2.1.1. Choice of modulus}
Our current goal is to find good values for the parameters that define
a linear congruential sequence.
Let us first consider the proper choice of the number $m$. We
want $m$ to be rather large, since the period cannot have more than
$m$ elements.\xskip (Even if a person wants
to generate only random zeros and ones, he should {\sl not}
take $m = 2$, for the sequence at best would then have the form
$\ldotss, 0, 1, 0, 1, 0, 1, \ldotss$! Methods for modifying random
numbers to get random zeros and ones are discussed in Section
3.4.)

Another factor that influences our choice of $m$ is speed of
generation: We want to pick a value so that the computation
of $(aX↓n + c)\mod m$ is fast.

Consider \MIX\ as an example. We can compute $y \mod m$ by putting
$y$ in registers A and X and dividing by $m$; assuming that
$y$ and $m$ are positive, we see that $y \mod m$ will then
appear in register X. But division is a comparatively slow operation,
and it can be avoided if we take $m$ to be a value that is
especially convenient, such as the {\sl word size} of our computer.

Let $w$ be the computer's word size, namely, $2↑e$ on an $e$-bit
binary computer or $10↑e$ on an $e$-digit decimal machine.\xskip (In
this book we shall often use the letter $e$ to denote an arbitrary
integer exponent, instead of the base of natural logarithms,
hoping that the context will make our notation unambiguous.
Physicists have a similar problem when they use $e$ for the
charge on an electron.)\xskip The result of an addition operation
is usually given modulo $w$, except on ones'-complement machines;
and multiplication mod $w$ is also quite simple, since the desired
result is the lower half of the product. Thus, the following
program computes the quantity $(aX + c)\mod w$ efficiently:
$$\vcenter{\mixtwo{
LDA⊗A⊗$\rA ←a$.\cr
MUL⊗X⊗$\rAX ← (\rA) \cdot X$.\cr
SLAX⊗5⊗$\rA ← \rAX \mod w$.\cr
ADD⊗C⊗$\rA ← (\rA + c)\mod w$.\quad\blackslug\cr}}\eqno (1)$$
The result appears in register A. The overflow toggle
might be on at the conclusion of the above sequence of instructions,
and if this is undesirable, the code should be followed by,
e.g., ``\.{JOV *+1}'' to turn it off.
%folio 13 galley 5 (C) Addison-Wesley 1978	*
A clever technique that is less commonly known can be used
to perform computations modulo $(w + 1)$. For reasons to be explained
later, we will generally want $c = 0$ when $m = w + 1$, so we
merely need to compute $(aX)\mod(w + 1)$. The following program
does this:
$$\vcenter{\mixfour{
01⊗⊗LDAN⊗X⊗$\rA←-X$.\cr
02⊗⊗MUL⊗A⊗$\rAX ← (\rA) \cdot a$.\cr
03⊗⊗STX⊗TEMP\cr
04⊗⊗SUB⊗TEMP⊗$\rA ← \rA - \rX$.\cr
05⊗⊗JANN⊗*+3⊗Exit if $\rA ≥ 0$.\cr
06⊗⊗INCA⊗2⊗$\rA←\rA + 2$.\cr
07⊗⊗ADD⊗=$w-1$=⊗$\rA←\rA + w - 1$.\xskip(Cf.\ exercise 3.)\quad\blackslug
\cr}}\eqno(2)$$
Register A now contains the value $(aX)\mod(w+1).$
Of course, this value might lie anywhere between 0 and
$w$, inclusive, so the reader may legitimately wonder how we
can represent so many values in the A-register! (The register
obviously cannot hold a number larger than $w - 1$.) The answer
is that overflow will be on after the above program if and only
if the result equals $w$, assuming that overflow was initially
off. It is convenient simply to reject the value $w$ if it
appears in the congruential sequence modulo $w + 1$; this will
happen if lines 05 and 06 of (2) are replaced by ``\.{JANN *+4;}
\.{INCA 2;} \.{JAP *-5}''.

To prove that code (2) actually does determine $(aX)\mod(w+1)$,
note that in line 04 we are subtracting the lower half
of the product from the upper half. No overflow can occur at
this step; and if $aX = qw + r$, with $0 ≤ r < w$, we will have
the quantity $r - q$ in register A after line 04. Now
$$aX = q(w + 1) + (r - q),$$ and since $q < w$,
we have $-w < r - q < w$; hence $(aX)\mod(w + 1)$ equals
either $r - q$ or $r - q + (w + 1)$, depending on whether $r
- q ≥ 0$ or $r - q < 0$.

A similar technique can be used to get the product
of two numbers modulo $(w - 1)$; see exercise 8.

In later sections we shall require a knowledge of the prime
factors of $m$ in order to choose the multiplier $a$ correctly.
Table 1 lists the complete factorization of $w \pm 1$ into primes
for nearly every known computer word size; the methods of Section 4.5.4 can be
used to extend this table if desired.

\topinsert{\vbox to 520pt{\hbox{(Table 1 will go on this page, 
it's being set separately)}}}

The reader may well ask why we bother to consider using $m =
w \pm 1$, when the choice $m = w$ is so manifestly convenient.
The reason is that {\sl when $m = w$, the right-hand digits of
$X↓n$ are much less random than the left-hand digits.} If $d$
is a divisor of $m$, and if
$$\eqalignno{Y↓n ⊗= X↓n \mod d,⊗(3)\cr
\noalign{\vskip 6pt\hbox{we can easily show that}\vskip 6pt}
Y↓{n+1} ⊗= (aY↓n + c)\mod d.⊗(4)\cr}$$
(For, $X↓{n+1} = aX↓n + c - qm$ for some integer $q$,
and taking both sides mod $d$ causes the quantity $qm$ to drop
out when $d$ is a factor of $m$.)

To illustrate the significance of Eq.\ (4), let
us suppose, for example, that we have a binary computer. If
$m = w = 2↑e$, the low-order four bits of $X↓n$ are the numbers
$Y↓n = X↓n \mod 2↑4$. The gist of Eq.\ (4) is that the low-order
four bits of $\langle X↓n \rangle$ form a congruential sequence that has a
period of length 16 or less.
Similarly, the low-order
five bits are periodic with a period of at most 32; and the
least significant bit of $X↓n$ is either constant or strictly
alternating.

This situation does not occur when $m = w \pm 1;$ in such a case,
the low-order bits of $X↓n$ will behave just as randomly as
the high-order bits do. If, for example, $w = 2↑{35}$ and $m
= 2↑{35} - 1$, the numbers of the sequence will not be very random
if we consider only their remainders mod 31, 71, 127, or 122921
(cf.\ Table 1); but the low-order bit, which represents the numbers
of the sequence taken mod 2, would be satisfactorily random.

Another alternative is to let $m$ be the largest prime number
less than $w$. This prime may be found by using the techniques
of Section 4.5.4, and a table of suitably large primes appears
in that section.

In most applications, the low-order bits are insignificant,
and the choice $m = w$ is quite satisfactory---provided that
the programmer using the random numbers does so wisely.
%folio 15 galley 6 (C) Addison-Wesley 1978	*
\yskip Our discussion so far has been based on
a ``signed magnitude'' computer like \MIX. Similar ideas apply
to machines that use complement notations, although there are
some instructive variations. For example, the PDP-10 computer
has 36 bits with two's complement arithmetic; when it computes
the product of two nonnegative integers, the lower half contains
the least significant 35 bits with a plus sign. On this machine
we should therefore take $w = 2↑{35}$, not $2↑{36}$. The 32-bit
two's complement arithmetic on IBM System/370 computers is different:
the lower half of a product contains a full 32 bits. Some programmers
have felt that this is a disadvantage, since the lower half
can be negative when the operands are positive, and it is a
nuisance to correct this; but actually it is a distinct {\sl
advantage} from the standpoint of random-number generation,
since we can take $m = 2↑{32}$ instead of $2↑{31}$ (see exercise
4).

\exbegin{EXERCISES}
\exno 1. [M12] In exercise
3.2.1--3 we concluded that the best congruential generators
will have $a$ relatively prime to $m$. Show that when $m = w$
in this case it is possible to compute $(aX + c)\mod w$ in
just {\sl three\/} \MIX\ instructions, rather than the four in (1),
with the result appearing in register X.

\exno 2. [16] Write a \MIX\ subroutine having the following characteristics:
$$\eqalign{\hbox{Calling sequence:\qquad}⊗\.{JMP RANDM}\cr
\noalign{\vskip 2pt}
\hbox{Entry conditions:\qquad}⊗
\hbox{Location \.{XRAND} contains an integer $X$.}\cr
\noalign{\vskip 2pt}
\hbox{Exit conditions:\qquad}⊗
\hbox{$X ← \rA ← (aX + c)\mod w$, $\rX ← 0$, overflow off.}\cr}$$
(Thus a call on this subroutine will produce the next
random number of a linear congruential sequence.)

\exno 3. [20] How can the constant $(w - 1)$ be specified in
general in the \MIX\ assembly language, regardless of the value
of the byte size?

\trexno 4. [21] Discuss the calculation of linear congruential
sequences with $m = 2↑{32}$ on two's-complement machines such
as the System/370 series.

\exno 5. [20] Given that $m$ is less than the word size, and
that $x, y$ are nonnegative integers less than $m$, show that
the difference $(x - y)\mod m$ may be computed in just four
\MIX\ instructions, without requiring any division. What is the
best code for the sum $(x + y)\mod m$?

\trexno 6. [20] The previous exercise
suggests that subtraction mod $m$ is easier to perform than
addition mod $m$. Discuss sequences generated by the rule
$$X↓{n+1} = (aX↓n - c)\mod m.$$ Are these sequences
essentially different from linear congruential sequences as
defined in the text? Are they more suited to efficient computer
calculation?

\exno 7. [M24] What patterns can you spot in Table 1?

\trexno 8. [20] Write a \MIX\ program analogous to (2) that computes
$(aX)\mod(w - 1)$. The values 0 and $w - 1$ are to be treated
as equivalent in the input and output of your program.

\exno 9. [23] Write a \MIX\ program analogous to the one in exercise
8, but it should compute $(aX)\mod(w - 2)$.

\runningrighthead{CHOICE OF MULTIPLIER} \section{3.2.1.2}
\dimsectionbegin{3.2.1.2. Choice of multiplier}
In this section we shall show how to choose the multiplier $a$
so as to give the {\sl period of maximum length.} A long period
is essential for any sequence that is to be used as a source
of random numbers; indeed, we would hope that the period contains
considerably more numbers than will ever be used in a single
application. Therefore we shall concern ourselves in this section
with the question of period length. The reader should keep in
mind, however, that a long period is only one desirable criterion
for the randomness of our sequence. For example, when $a = c
= 1$, the sequence is simply $X↓{n+1} = (X↓n + 1)\mod m$, and
this obviously has a period of length $m$, yet it is anything
but random. Other considerations affecting the choice of a multiplier
will be given later in this chapter.

Since only $m$ different values are possible, the period cannot
be longer than $m$. Can we achieve the maximum length, $m$?
The example above shows that it is always possible, although
the choice $a = c = 1$ does not yield a desirable sequence.
Let us investigate {\sl all} possible choices of $a$, $c$, and
$X↓0$ that give a period of length $m$. It turns out that all
such values of the parameters can be characterized very simply;
when $m$ is the product of distinct primes, only $a = 1$ will
produce the full period, but when $m$ is divisible by a high
power of some prime there is considerable latitude in the choice
of $a$.

\thbegin Theorem A. {\sl The linear congruential
sequence has period length $m$ if and only if
$$\baselineskip 15pt
\vbox{\halign to size{\rt{\qquad\rm #\xskip}⊗\lft{#}\tabskip0pt plus 300pt\cr
i)⊗$c$ is relatively prime to $m$;\cr
ii)⊗$b = a - 1$ is a multiple of $p$, for every prime $p$ dividing $m$;\cr
iii)⊗$b$ is a multiple of $4$, if $m$ is a multiple of $4$.\cr}}$$}	%end \sl

The ideas used in the proof of this theorem
go back at least a hundred years. The first proof of the theorem
in this particular form was given by M. Greenberger in the special
case $m = 2↑e$ [see {\sl JACM \bf 8} (1961), 383--389],
and the sufficiency of conditions (i), (ii), and (iii) in
the general case was shown by Hull and Dobell [see {\sl
SIAM Review \bf 4} (1962), 230--254]. To prove the theorem
we will first consider some auxiliary number-theoretic results
that are of interest in themselves.

\thbegin Lemma P. {\sl Let $p$ be a prime number,
and let $e$ be a positive integer, where $p↑e > 2$. If}
$$\eqalignno{x ≡ 1 \modulo{p↑e},⊗\qquad x \neqv 1 \modulo
{p↑{e+1}}⊗(1)\cr\noalign{\hbox{\sl then}}
x↑p ≡ 1 \modulo{p↑{e+1}},⊗\qquad x↑p \neqv 1 \modulo
{p↑{e+2}}.⊗(2)\cr}$$

\dproofbegin We
have $x = 1 + qp↑e$ for some integer $q$ that is not a multiple
of $p$. By the binomial formula
$$\eqalign{x↑p⊗=1+{p\choose 1}qp↑e+\cdots+{p\choose{p-1}}q↑{p-1}
p↑{(p-1)e}+q↑p p↑{pe}\cr
⊗=1+qp↑{e+1}\left(1+{1\over p}{p\choose 2}qp↑e + {1\over p}
{p\choose 3}q↑2 p↑{2e}+\cdots+{1\over p}{p\choose p}q↑{p-1}
p↑{(p-1)e}\right).\cr}$$
The quantity in parentheses is an integer, and, in fact,
every term in the parentheses is a multiple of $p$ except the
first term. For if $1 < k < p$, the binomial coefficient $p \choose k$
is divisible by $p$ (cf.\ exercise 1.2.6--10), hence
$${1\over p}{p \choose k} q↑{k-1} p↑{(k-1)e}$$
is divisible by $p↑{(k-1)e}$; and the last term is
$q↑{p-1}p↑{(p-1)e-1}$, which is divisible by $p$ since
$(p - 1)e > 1$ when $p↑e > 2$. So $x↑p ≡ 1 + qp↑{e+1} \modulo
{p↑{e+2}}$, and this completes the proof.\xskip $\biglp${\sl Note:}
A generalization of this result appears in exercise 3.2.2--11(a).$\bigrp$
\quad \blackslug

\thbegin Lemma Q. {\sl Let the decomposition
of $m$ into prime factors be
$$m = p↑{e↓1}↓{1} \ldotss p↑{e↓t}↓{t}.\eqno (3)$$
The length, $λ$, of the period of the linear congruential sequence
determined by $(X↓0, a, c, m)$ is the least common multiple of
the lengths $λ↓j$ of the periods of the linear congruential sequences
$(X↓0\mod p↑{e↓j}↓j, a\mod p↑{e↓j}↓j, c\mod p↑{e↓j}↓j,
p↑{e↓j}↓{j})$, $1 ≤ j ≤ t$.}
%folio 20 galley 7 (C) Addison-Wesley 1978	*
\proofbegin By induction on $t$, it suffices
to prove that if $m↓1$ and $m↓2$ are relatively prime, the length
$λ$ of the linear congruential sequence determined by $(X↓0,
a, c, m↓1m↓2)$ is the least common multiple of the lengths $λ↓1,
λ↓2$ of the periods of the sequences $(X↓0\mod m↓1,\, a\mod m↓1,\,
c\mod m↓1,\,m↓1)$ and $(X↓0\mod m↓2,\,a\mod m↓2,\, c\mod m↓2,\,m↓2)$.
We observed in the previous section, Eq.\ (4), that if the elements
of these three sequences are denoted by $X↓n$, $Y↓n$, and $Z↓n$,
respectively, we will have
$$Y↓n = X↓n\mod m↓1\qquad\hbox{and}\qquad Z↓n = X↓n\mod m↓2
\qquad\hbox{for all }n ≥ 0.$$
Therefore, by Law D of Section 1.2.4, we find that
$$X↓n = X↓k\qquad\hbox{if and only if}\qquad Y↓n = Y↓k\quad
\hbox{and}\quad Z↓n = Z↓k.\eqno (4)$$

Let $λ↑\prime$
be the least common multiple of $λ↓1$ and $λ↓2$; we wish to
prove that $λ↑\prime = λ$. Since $X↓n = X↓{n+λ}$ for all suitably
large $n$, we have $Y↓n = Y↓{n+λ}$ (hence $λ$ is a multiple
of $λ↓1$) and $Z↓n = Z↓{n+λ}$ (hence $λ$ is a multiple of $λ↓2$),
so $λ$ must be $≥λ↑\prime$. Furthermore, we know that $Y↓n
= Y↓{n+λ↑\prime}$ and $Z↓n = Z↓{n+λ↑\prime}$ for all suitably
large $n;$ therefore, by (4), $X↓n = X↓{n+λ↑\prime}$. This
proves $λ ≤ λ↑\prime$, so $λ = λ↑\prime$.\quad\blackslug

\yyskip Now we are ready to prove Theorem A\null. Because
of Lemma Q\null, it suffices to prove the theorem when $m$ is a power
of a prime number. For
$$p↑{e↓1}↓{1} \ldotss p↑{e↓t}↓{t} = λ = \hbox{lcm}(λ↓1, \ldotss, λ↓t)
≤ λ↓1 \ldotsm λ↓t ≤ p↑{e↓1}↓{1} \ldotss p↑{e↓t}↓{t}$$
can be true if and only if $λ↓j = p↑{e↓j}↓{\!j}$ for $1 ≤ j ≤ t$.

Therefore, assume that $m = p↑e$, where $p$ is
prime and $e$ is a positive integer. The theorem is obviously
true when $a = 1$, so we may take $a > 1$. The period can be
of length $m$ if and only if each possible integer $0 ≤ x < m$
occurs in the period, since no value occurs in the period
more than once. Therefore the period is of length $m$ if and
only if the period of the sequence with $X↓0 = 0$ is of length
$m$, and we are justified in supposing that $X↓0 = 0$. By formula
3.2.1--6 we have
$$X↓n = \left({a↑n - 1\over a - 1}\right)\,c\, \mod\, m.\eqno (5)$$
If $c$ is not relatively prime to $m$, this value $X↓n$ could never
be equal to 1, so condition (i) of the theorem is necessary. The period has
length $m$ if and only if the smallest positive value of $n$
for which $X↓n = X↓0 = 0$ is $n = m$. By (5) and condition (i),
our theorem now reduces to proving the following fact:

\thbegin Lemma R. {\sl Assume that $1 < a < p↑e$, where
$p$ is prime. If $λ$ is the smallest positive integer for which
$(a↑λ - 1)/(a - 1) ≡ 0 \modulo {p↑e}$, then
$$λ = p↑e\qquad\hbox{if and only if}\qquad \left\{
\vcenter{\halign{\lft{$#$}⊗\qquad#\quad⊗\lft{$#$}\cr
a ≡ 1 \modulo{p}⊗when⊗p>2,\cr
a ≡ 1 \modulo{4}⊗when⊗p=2.\cr}}\right.$$}

\dproofbegin Assume that $λ = p↑e$. If $a \neqv 1 \modulo {p}$, then
$(a↑n - 1)/(a - 1) ≡ 0 \modulo {p↑e}$ if and only if $a↑n -
1 ≡ 0 \modulo {p↑e}$. The condition $a↑{p↑e} - 1 ≡ 0
\modulo {p↑e}$ then implies that $a↑{p↑e} ≡ 1 \modulo{p}$;
but by Theorem 1.2.4F\null, $a↑{p↑e} ≡ a \modulo {p},$
hence $a \neqv 1 \modulo{p}$ leads to a contradiction.
And if $p = 2$ and $a ≡ 3 \modulo{4}$, we have $(a↑{2
↑{e-1}} - 1)/(a - 1) ≡ 0 \modulo {2↑e}$ by exercise
8. These arguments show that it is necessary in general to have $a = 1
+ qp↑f$, where $p↑f > 2$ and $q$ is not a multiple of $p$, whenever
$λ = p↑e.$

It remains to be shown that this condition is {\sl sufficient}
to make $λ = p↑e$. By repeated application of Lemma P\null, we find
that
$$a↑{p↑g} ≡ 1 \modulo {p↑{f+g}},\qquad a↑{p↑g} \neqv 1 \modulo {p↑{f+g+1}},$$
for all $g ≥ 0$, and therefore
$$\eqalign{(a↑{p↑g} - 1)/(a - 1)⊗≡ 0 \modulo {p↑g},\cr
(a↑{p↑g} - 1)/(a - 1)⊗\neqv 0 \modulo {p↑{g+1}}.\cr}\eqno(6)$$
In particular, $(a↑{p↑e} - 1)/(a - 1)
≡ 0 \modulo {p↑e}$. Now the congruential sequence $(0, a, 1,
p↑e)$ has $X↓n = (a↑n - 1)/(a - 1)\mod p↑e$; therefore it
has a period of length $λ$, that is, $X↓n = 0$ if and only if
$n$ is a multiple of $λ$. Hence $p↑e$ is a multiple of $λ$.
This can happen only if $λ = p↑g$ for some $g$, and the relations
in (6) imply that $λ = p↑e$, completing the proof.\quad\blackslug

The proof of Theorem A is now complete.\quad\blackslug

\yskip We will conclude this section by considering
the special case of pure multiplicative generators, when $c
= 0$. Although the random-number generation process is slightly
faster in this case, Theorem A shows us that the maximum period
length cannot be achieved. In fact, this is quite obvious, since
the sequence now satisfies the relation
$$X↓{n+1} = aX↓n\mod m,\eqno (7)$$ and the value
$X↓n = 0$ should never appear lest the sequence degenerate to
zero. In general, if $d$ is any divisor of $m$ and if $X↓n$
is a multiple of $d$, all succeeding elements $X↓{n+1}$, $X↓{n+2}$,
$\ldots$ of the multiplicative sequence will be multiples of
$d$. So when $c = 0$, we will want $X↓n$ to be relatively prime
to $m$ for all $n$, and this limits the length of the period to
at most $\varphi(m)$, the number of integers between 0 and $m$ that
are relatively prime to $m$.

It may be possible to achieve an acceptably long period even
if we stipulate that $c = 0$. Let us now try to find conditions
on the multiplier so that the period is as long as possible
in this special case.

According to Lemma Q\null, the period of the sequence depends entirely
on the periods of the sequences when $m = p↑e$, so let us consider
that situation. We have $X↓n = a↑nX↓0\mod p↑e$, and it is
clear that the period will be of length 1 if $a$ is a multiple
of $p$, so we take $a$ to be relatively prime to $p$. Then
the period is the smallest integer $λ$ such that $X↓0 = a↑λX↓0
\mod p↑e$. If the greatest common divisor of $X↓0$ and $p↑e$
is $p↑f$, this condition is equivalent to
$$a↑λ ≡ 1 \modulo {p↑{e-f}}.\eqno (8)$$
By Euler's theorem (exercise 1.2.4--28), $a↑{\varphi(p↑{e-f})}
≡ 1 \modulo {p↑{e-f}};$ hence $λ$ is a divisor of
$$\varphi(p↑{e-f}) = p↑{e-f-1}(p - 1).$$
When $a$ is relatively prime to $m$, the smallest integer $λ$
for which $a↑λ ≡ 1 \modulo {m}$ is conventionally called {\sl
the order of $a$ modulo $m$.} Any such value of $a$ that gives
the {\sl maximum} possible order modulo $m$ is called a {\sl
primitive element} modulo $m$.

Let $λ(m)$ denote the order of a primitive element, i.e., the
maximum possible order, modulo $m$. The remarks above show that
$λ(p↑e)$ is a divisor of $p↑{e-1}(p - 1);$ with a little care
(see exercises 11 through 16 below) we can give the precise
value of $λ(m)$ in all cases as follows:
$$\vcenter{\halign{\ctr{#}\cr
$λ(2) = 1$,\qquad $λ(4) = 2$,\qquad $λ(2↑e) = 2↑{e-2}$\quad if\quad $e≥3$.\cr
$λ(p↑e) = p↑{e-1}(p - 1)$,\qquad if\quad $p > 2$.\cr
$λ(p↑{e↓1}↓{1} \ldotss p↑{e↓t}↓{t}) = \hbox{lcm}\left( λ(p↑{e↓1}↓{1}),
\ldotss, λ(p↑{e↓t}↓{t})\right)$.\cr}}\eqno(9)$$ Our remarks may
be summarized in the following theorem:

\thbegin Theorem B. [R. D. Carmichael, {\sl Bull.\ Amer.\ 
Math.\ Soc.\ \bf 16} (1910), 232--238.] {\sl The maximum period possible
when $c = 0$ is $λ(m)$, where $λ(m)$ is defined in $(9)$. This period
is achieved if
$$\baselineskip 15pt
\vbox{\halign to size{\rt{\qquad\rm#\xskip}⊗\lft{#}\tabskip0pt plus 300pt\cr
i)⊗$X↓0$ is relatively prime to $m$;\cr
ii)⊗$a$ is a primitive element modulo $m$.\quad\blackslug\cr}}$$}	%end \sl

\noindent Note that we can obtain a period of length $m-1$
if $m$ is prime; this is just one less than the maximum length, so for
all practical purposes the period is as long as we want.

The question now is, how can we find primitive elements modulo
$m$? The exercises at the close of this section give us the
following facts:
%folio 23 galley 8 (C) Addison-Wesley 1978	*
\thbegin Theorem C. {\sl The number $a$ is a primitive
element modulo $p↑e$ if and only if
$$\vbox{\halign to size{\rt{\qquad\rm#\xskip}⊗\lft{#}\tabskip0pt plus 300pt\cr
i)⊗$p↑e = 2$, $a$ is odd;\qquad or $p↑e = 4$, $a\mod 4 = 3$;\cr
⊗or $p↑e = 8$, $a\mod 8 = 3,5,7$;\qquad or
$p=2$, $e≥4$, $a\mod 8 = 3, 5$;\cr
\noalign{\yskip\hbox{or}\yskip}
ii)⊗$p$ is odd, $e=1$, $a \neqv 0 \modulo{p}$, and
$a↑{(p-1)/q}\neqv 1 \modulo{p}$\cr
⊗for any prime divisor $q$ of $p-1$;\cr
\noalign{\yskip\hbox{or}\yskip}
iii)⊗$p$ is odd, $e>1$, $a$ satisfies (ii),
and $a↑{p-1} \neqv 1 \modulo{p↑2}$.\quad\blackslug\cr}}$$}	%end \sl

\noindent Conditions (ii) and (iii) of this
theorem are readily tested on a computer for large values of
$p$, by using the efficient methods for evaluating powers discussed
in Section 4.6.3. Theorem C applies to powers of primes only;
if we are given values $a↓j$ that
are primitive modulo $p↑{e↓j}↓{\!j}$, it is possible to find a
single value $a$ such that $a ≡ a↓j \modulo {p↑{e↓j}↓{\!j}}$, for
$1 ≤ j ≤ t$, using the ``Chinese remainder algorithm'' discussed
in Section 4.3.2, and this number $a$ will be a primitive element modulo
$p↑{e↓1}↓{1} \ldotss p↑{e↓t}↓{t}$. Hence there is a reasonably
efficient way to construct multipliers satisfying the condition
of Theorem B\null, for any desired value of $m$, although the calculations
can be somewhat lengthy in the general case.

In the common case $m = 2↑e$, with $e ≥ 4$, the conditions above
simplify to the single requirement that $a ≡ 3$ or $5 \modulo
{8}$. In this case, one-fourth of all possible multipliers give
the maximum period.

The second most common case is when $m = 10↑e$. Using Lemmas
P and Q\null, it is not difficult to obtain necessary and sufficient
conditions for the achievement of the maximum period in the
case of a decimal computer (cf.\ exercise 18):

\thbegin Theorem D. {\sl If $m = 10↑e$, $e ≥ 5$, $c = 0$, and
$X↓0$ is not a multiple of $2$ or $5$, the period of the linear congruential
sequence is $5 \times 10↑{e-2}$ if and only if $a\mod 200$ equals
one of the following $32$ values:}
$$\vcenter{\halign{\quad\lft{#}\cr
3, 11, 13, 19,
21, 27, 29, 37, 53, 59, 61, 67, 69, 77, 83, 91, 109, 117,\cr
123, 131, 133, 139, 141, 147, 163, 171, 173, 179, 181,
187, 189, 197.\quad\blackslug\cr}}\eqno(10)$$

\exbegin{EXERCISES}
\exno 1. [10] What is
the length of the period of the linear congruential sequence
with $X↓0 = 5772156648$, $a = 3141592621$, $c = 2718281829$, and
$m = 10000000000$?

\exno 2. [10] Are the following two conditions sufficient to
guarantee the maximum length period, when $m$ is a power
of 2? ``(i) $c$ is odd; (ii) $a\mod 4 = 1$.''

\exno 3. [13] Suppose that $m = 10↑e$, where $e ≥ 2$, and suppose
further that $c$ is odd and not a multiple of 5. Show
that the linear congruential sequence will have the maximum
length period if and only if $a\mod 20 = 1$.

\exno 4. [M20] When $a$ and $c$ satisfy the conditions of Theorem
A\null, and when $m = 2↑e$, $X↓0 = 0$, what is the value of $X↓{2↑{e-1}}$?

\exno 5. [14] Find all multipliers $a$ that satisfy the conditions
of Theorem A when $m = 2↑{35} + 1$.\xskip (The prime factors of $m$
may be found in Table 3.2.1.1--1.)

\trexno 6. [20] Find all multipliers $a$ that satisfy the conditions
of Theorem A when $m = 10↑6 - 1$.\xskip (See Table 3.2.1.1--1.)

\trexno 7. [M23] The period of a congruential sequence need not
start with $X↓0$, but we can always find indices $\mu ≥ 0$ and $λ
> 0$ such that $X↓{n+λ} = X↓n$ whenever $n ≥ \mu $, and for
which $\mu$ and $λ$ are the smallest possible values with this
property.\xskip (Cf.\ exercises 3.1--6 and 3.2.1--1.)\xskip If $\mu ↓j, λ↓j$
are the indices corresponding to the sequences $(X↓0\mod p↑{e↓j}↓{\!j},\penalty0
a\mod p↑{e↓j}↓{\!j}, c\mod p↑{e↓j}↓{\!j}, p↑{e↓j}↓{\!j})$, and
if $\mu , λ$ correspond to the sequence $(X↓0, a, c, p↑{e↓1}↓{1}
\ldotss p↑{e↓t}↓{t})$, Lemma Q states that $λ$ is the least common
multiple of $λ↓1, \ldotss, λ↓t$. What is the value of $\mu$ in
terms of the values of $\mu ↓1, \ldotss, \mu ↓t$? What is the
maximum possible value of $\mu$ obtainable by varying $X↓0$,
$a$, and $c$, when $m = p↑{e↓1}↓{1} \ldotss p↑{e↓t}↓{t}$ is fixed?

\exno 8. [M20] Show that if $a\mod 4 = 3$, we have $(a↑{2↑{e-1}}
 - 1)/(a - 1) ≡ 0 \modulo {2↑e}$ when $e > 1$.\xskip (Use
Lemma P.)

\trexno 9. [M22] (W. E. Thomson.)\xskip When $c = 0$ and $m = 2↑e ≥
16$, Theorems B and C say that the period has length $2↑{e-2}$
if and only if the multiplier $a$ satisfies $a\mod 8 = 3$ or
$a\mod 8 = 5$. Show that every such sequence is essentially
a linear congruential sequence with $m = 2↑{e-2}$, having {\sl
full} period, in the following sense:

\yskip\textindent{a)}If $X↓{n+1} = (4c + 1)X↓n\mod 2↑e$, and
$X↓n = 4Y↓n + 1$, then$$
Y↓{n+1} = \biglp (4c + 1)Y↓n + c\bigrp\mod
2↑{e-2}.$$

\textindent{b)}If $X↓{n+1}
= (4c - 1)X↓n \mod 2↑e$, and $X↓n = \biglp (-1)↑n(4Y↓n +
1)\bigrp\mod 2↑e$, then
$$Y↓{n+1} = \biglp (1 - 4c)Y↓n - c\bigrp\mod 2↑{e-2}.$$

({\sl Note:} In these formulas, $c$ is
an odd integer. The literature contains several statements to
the effect that sequences with $c = 0$ satisfying Theorem B
are somehow more random than sequences satisfying Theorem A\null,
in spite of the fact that the period is only one-fourth as long
in the case of Theorem B\null. This exercise refutes such statements;
in essence, one gives up two bits of the word length in order to
save the addition of $c$, when $m$ is a power of 2.)

\exno 10. [M21] For what values of $m$ is $λ(m) = \varphi
(m)$?

\trexno 11. [M28] Let $x$ be an odd integer greater than 1.\xskip (a)
Show that there exists a unique integer $f > 1$ such that $x
≡ 2↑f \pm 1 \modulo{2↑{f+1}}.$\xskip (b) Given that $1 < x < 2↑e
- 1$ and that $f$ is the corresponding integer from part (a),
show that the order of $x$ modulo $2↑e$ is $2↑{e-f}$.\xskip (c) In
particular, this proves Theorem C(i).

\exno 12. [M26] Let $p$ be an odd prime. If $e > 1$, prove that
$a$ is a primitive element modulo $p↑e$ if and only if $a$ is
a primitive element modulo $p$ and $a↑{p-1} \neqv 1 \modulo
{p↑2}$.\xskip $\biglp$For the purposes of this exercise, assume that
$λ(p↑e) = p↑{e-1}(p - 1)$. This fact is proved in exercises
14 and 16 below.$\bigrp$

\exno 13. [M22] Let $p$
be prime. Given that $a$ is not a primitive element modulo $p$,
show that either $a$ is a multiple of $p$ or $a↑{(p-1)/q}
≡ 1 \modulo {p}$ for some prime number $q$ that divides $p-1$.

\exno 14. [M18] If $e > 1$ and $p$ is an odd prime, and if $a$
is a primitive element modulo $p$, prove that either $a$ or
$a + p$ is a primitive element modulo $p↑e$.\xskip [{\sl Hint:} See
exercise 12.]

\exno 15. [M29] (a) Let $a↓1, a↓2$ be relatively prime to $m$,
and let their orders modulo $m$ be $λ↓1, λ↓2$, respectively.
If $λ$ is the least common multiple of $λ↓1$ and $λ↓2$, prove
that $a↑{\kappa ↓1}↓{1} a↑{\kappa ↓2}↓{2}$ has order $λ$ modulo
$m$, for suitable integers $\kappa ↓1, \kappa ↓2$.\xskip [{\sl Hint:}
Consider first the case that $λ↓1$ is relatively prime to
$λ↓2.$]\xskip (b) Let $λ(m)$ be the maximum order of any element modulo
$m$. Prove that $λ(m)$ is a multiple of the order of each element
modulo $m$; that is, prove that $a↑{λ(m)} ≡ 1 \modulo {m}$
whenever $a$ is relatively prime to $m$.

\trexno 16. [M24] Let $p$ be a prime number.\xskip (a) Let $f(x) =
x↑n + c↓1x↑{n-1} + \cdots + c↓n$, where the $c$'s
are integers. Given that $a$ is an integer for which $f(a) ≡
0 \modulo {p}$, show that there exists a polynomial $q(x) =
x↑{n-1} + q↓1x↑{n-2} + \cdots + q↓{n-1}$ with integer
coefficients such that $f(x) ≡ (x - a)q(x) \modulo {p}$ for
all integers $x$.\xskip (b) Let $f(x)$ be a polynomial as in (a).
Show that $f(x)$ has at most $n$ distinct ``roots'' modulo $p$;
that is, there are at most $n$ integers $a$, with $0 ≤ a < p$,
such that $f(a) ≡ 0 \modulo {p}$.\xskip (c) Because
of exercise 15(b), the polynomial $f(x) = x↑{λ(p)} - 1$ has
$p - 1$ distinct roots; hence there is an integer $a$ with order
$p - 1$.

\exno 17. [M26] Not all of the values listed in Theorem D would
be found by the text's construction; for example, 11 is not
primitive modulo $5↑e$. How can this be possible, when 11 {\sl is}
primitive modulo $10↑e$, according to Theorem D\null?
Which of the values listed in Theorem D are primitive
elements modulo {\sl both} $2↑e$ and $5↑e$?

\exno 18. [M25] Prove Theorem D\null.\xskip (Cf.\ the previous exercise.)

\exno 19. [40] Make a table of some suitable multipliers, $a$,
for each of the values of $m$ listed in Table 3.2.1.1--1, assuming
that $c = 0$.

\trexno 20. [M24] (G. Marsaglia.)\xskip The purpose of this exercise
is to study the period length of an {\sl arbitrary} linear congruential
sequence. Let $Y↓n = 1 + a + \cdots + a↑{n-1}$, so
that $X↓n = (AY↓n + X↓0)\mod m$ for some constant $A$ by Eq.\ 3.2.1--8.\xskip
(a) Prove that the period length of $\langle X↓n\rangle$ is
the period length of $\langle Y↓n \mod m↑\prime \rangle $, where
$m↑\prime = m/\!\gcd(A, m)$.\xskip (b) Prove that the period length
of $\langle Y↓n \mod p↑e\rangle$ satisfies the following when
$p$ is prime:\xskip (i) If $a \mod p = 0$, it is 1.\xskip (ii) If $a \mod
p = 1$, it is $p↑e$, except when $p = 2$ and $e ≥ 2$ and $a
\mod 4 = 3.$\xskip (iii) If $p = 2,\ e ≥ 2$, and $a \mod 4 = 3$, it is
twice the order of $a$ modulo $p↑e$ (cf.\ exercise 11), unless $a ≡ -1
\modulo{2↑e}$ when it is 2.\xskip (iv) If $a\mod p > 1$,
it is the order of $a$ modulo $p↑e$.

\exno 21. [M25] In a linear congruential sequence of maximum period, let
$X↓0 = 0$ and let $s$ be the least positive integer such that $a↑s ≡ 1
\modulo{m}$. Prove that $\gcd(X↓s,m) = s$.
%folio 26 galley 9 (C) Addison-Wesley 1978	*
\runningrighthead{POTENCY} \section{3.2.1.3}
\dimsectionbegin{3.2.1.3. Potency} In the preceding
section, we showed that the maximum period can be obtained when
$b = a - 1$ is a multiple of each prime dividing $m$; and $b$
must also be a multiple of 4 if $m$ is a multiple of 4. If
$z$ is the radix of the machine being used---so that $z = 2$
for a binary computer, and $z = 10$ for a decimal computer---and
if $m$ is the word size $z↑e$, the multiplier
$$a = z↑k + 1,\qquad 2 ≤ k < e\eqno (1)$$
satisfies these conditions. Theorem 3.2.1.2A also says that we may take
$c = 1$. The recurrence relation now has the form
$$X↓{n+1} = \biglp (z↑k + 1)X↓n + 1\bigrp \mod z↑e,\eqno(2)$$
and this equation suggests that we can avoid
the multiplication; merely shifting and adding will suffice.

For example, suppose that $a = B↑2 + 1$, where $B$ is the byte
size of \MIX. The code
$$\vcenter{\mixtwo{
LDA⊗X\cr
SLA⊗2\cr
ADD⊗X\cr
INCA⊗1\cr}}\eqno(3)$$
can be used in place of the instructions given in Section 3.2.1.1,
and the execution time decreases from $16u$ to $7u$.

For this reason, multipliers having form (1) have been widely
discussed in the literature, and indeed they have been recommended
by many authors. However, the early years of experimentation
with this method showed that {\sl multipliers having the simple
form in $(1)$ should be avoided.} The generated numbers just aren't
random enough.

Later in this chapter we shall be discussing some rather sophisticated
theory that accounts for the badness of all the linear congruential
random-number generators known to be bad. However, some generators
$\biglp$such as (2)$\bigrp$ are sufficiently awful that a comparatively
simple theory can be used to dispense with them. This simple
theory is related to the concept of ``potency,'' which we shall
now discuss.

The {\sl potency} of a linear congruential sequence with maximum
period is defined to be the least integer $s$ such that
$$b↑s ≡ 0 \modulo {m}.\eqno (4)$$
(Such an integer $s$ will always exist when the multiplier satisfies the
conditions of Theorem 3.2.1.2A\null, since $b$ is a multiple of every
prime dividing $m$.)

We may analyze the randomness of the sequence
by taking $X↓0 = 0$, since 0 occurs somewhere in the period.
With this assumption, we have $$X↓n = \biglp (a↑n - 1)c/b\bigrp
\mod m,$$ and if we expand $a↑n - 1 = (b + 1)↑n - 1$ by the
binomial theorem, we find that
$$X↓n = c\,\left(n + {n\choose 2}\,b + \cdots + {n\choose s}\,b↑{s-1}\right)
\mod m.\eqno (5)$$ All terms in $b↑s$,
$b↑{s+1}$, etc., may be ignored, since they are multiples of
$m$.

Equation (5) can be instructive, so we shall consider
some special cases. If $a = 1$, the potency is 1; and $X↓n ≡
cn \modulo {m}$, as we have already observed, so the sequence
is surely not random. If the potency is 2, we have $X↓n ≡ cn
+ cb{n\choose 2}$, and again the sequence is not very random;
indeed,
$$X↓{n+1} - X↓n ≡ c + cbn,$$
so the differences between consecutively generated numbers change in a very
simple way from one value of $n$ to the next. The point $(X↓n, X↓{n+1}, X↓{n+2})$
always lies on one of the four planes
$$\eqalign{x-2y+z⊗=d+m,\cr
x-2y+z⊗=d,\cr
x-2y+z⊗=d-m,\cr
x-2y+z⊗=d-2m,\cr}$$
in three-dimensional space, where
$d = cb \mod m$.

If the potency is 3, the sequence begins to look somewhat more
random, but there is a high degree of dependency between $X↓n$,
$X↓{n+1}$, and $X↓{n+2}$; tests show that sequences with potency 3 are still not
sufficiently good. Reasonable results have been reported when
the potency is 4 or more, but these have been disputed by other
people. A potency of at least 5 would seem to be required for
sufficiently random values.

Suppose, for example, that $m = 2↑{35}$ and $a = 2↑k + 1$. Then
$b = 2↑k$, so we find that when $k ≥ 18$, the value $b↑2 = 2↑{2k}$
is a multiple of $m$: the potency is 2. If $k = 17, 16, \ldotss,
12$, the potency is 3, and a potency of 4 is achieved for
$k = 11, 10, 9$. The only acceptable multipliers, from the standpoint
of potency, therefore have $k ≤ 8$. This means $a ≤ 257$, and
we shall see later that {\sl small} multipliers are also to
be avoided. We have now eliminated all multipliers of the
form $2↑k + 1$ when $m = 2↑{35}$.

When $m$ is equal to $w \pm 1$, where $w$ is the
word size, $m$ is generally not divisible by high powers of
primes, and a high potency is impossible (see exercise 6).
So in this case, the maximum-period method should {\sl not}
be used; the pure-multiplication method with $c = 0$ should
be applied instead.

It must be emphasized that high potency is necessary but not
sufficient for randomness; we use the concept of potency only
to reject impotent generators, not to accept the potent ones.
Linear congruential sequences should pass the ``spectral test''
of Section 3.3.4 before they are considered to be acceptably
random.

\exbegin{EXERCISES}
\exno 1. [M10] Show
that, no matter what the byte size $B$ of \MIX\ happens to be,
the code (3) yields a random-number generator of maximum period.

\exno 2. [10] What is the potency of the generator represented
by the \MIX\ code (3)?

\exno 3. [11] When $m = 2↑{35}$, what is the potency of the linear
congruential sequence with $a = 3141592621$? What is the potency
if the multiplier is $a = 2↑{23} + 2↑{14} + 2↑2 + 1$?

\exno 4. [15] Show that if $m
= 2↑e ≥ 8$, maximum potency is achieved when $a \mod 8 =
5.$

\exno 5. [M20] Given that $m = p↑{e↓1}↓{1} \ldotss p↑{e↓t}↓{t}$
and $a = 1 + kp↑{f↓1}↓{1} \ldotss p↑{f↓t}↓{t}$, where $a$ satisfies
the conditions of Theorem 3.2.1.2A and $k$ is relatively prime
to $m$, show that the potency is $\max(\lceil e↓1/f↓1\rceil ,
\ldotss, \lceil e↓t/f↓t\rceil )$.

\trexno 6. [20] Which of the values
of $m = w \pm 1$ in Table 3.2.1.1--1 can be used in a linear congruential
sequence of maximum period whose potency is 4 or more?
(Use the result of exercise 5.)

\exno 7. [M20] When $a$ satisfies the conditions of Theorem
3.2.1.2A\null, it is relatively prime to $m$; hence there is a number
$a↑\prime$ such that $aa↑\prime ≡ 1 \modulo {m}$. Show that $a↑\prime$
can be expressed simply in terms of $b$.

\trexno 8. [M26] A random-number
generator with $X↓{n+1} = (2↑{17} + 3)X↓n \mod 2↑{35}$ and $X↓0=1$
was subjected to the following test: Let $Y↓n = \lfloor 10X↓n/2↑{35}\rfloor
$; then $Y↓n$ should be a random digit between 0 and 9, and
the triples $(Y↓{3n},\ Y↓{3n+1},\ Y↓{3n+2})$ should take on each
of the 1000 possible values from (0, 0, 0) to (9, 9, 9) with
equal probability. But with 30000 values of $n$ tested, some
triples hardly ever occurred, and others occurred much more
often than they should have. Can you account for this failure?

\runningrighthead{OTHER METHODS} \section{3.2.2}
\sectionskip
\sectionbegin{3.2.2. Other Methods}
Of course, linear congruential sequences are
not the only sources of random numbers that have been proposed
for computer use. In this section we shall review the most significant
alternatives; some of these methods are quite important, while others
are interesting chiefly because they are not as good as a person
might expect.

One of the common fallacies encountered in connection with random-number
generation is the idea that we can take a good generator
and modify it a little, in order to get an ``even-more-random''
sequence. This is often false. For example, we know that
$$X↓{n+1} = (aX↓n + c)\mod m\eqno (1)$$ leads to
reasonably good random numbers; wouldn't the sequence produced by
$$X↓{n+1} = \biglp (aX↓n)\mod(m + 1) + c\bigrp \mod m\eqno(2)$$
be even {\sl more} random? The answer is, the
new sequence is probably a great deal {\sl less} random. For
the whole theory breaks down, and in the absence of any theory
about the behavior of the sequence (2), we come into the area
of generators of the type $X↓{n+1} = f(X↓n)$ with the function
$f$ chosen at random; exercises 3.1--11 through 3.1--15 show
that these sequences probably behave much more poorly than the sequences
obtained from the more disciplined function (1).

Let us consider another approach, in an attempt to get ``more
random'' numbers. The linear congruential method can be generalized
to, say, a quadratic congruential method:
$$X↓{n+1} = (dX↑2↓n + aX↓n + c)\mod m.\eqno (3)$$
Exercise 8 generalizes Theorem 3.2.1.2A to obtain necessary
and sufficient conditions on $a$, $c$, and $d$ such that the sequence
defined by (3) has a period of the maximum length $m$; the restrictions
are not much more severe than in the linear method.
%folio 32 galley 10 (C) Addison-Wesley 1978	*
An interesting quadratic method has been proposed by R. R. Coveyou
when $m$ is a power of two; let
$$X↓0\mod 4 = 2,\qquad X↓{n+1} = X↓n(X↓n + 1)\mod
2↑e,\qquad n ≥ 0.\eqno (4)$$
This sequence can be computed with about the same
efficiency as (1), without any worries of overflow. It has an
interesting connection with von Neumann's original middle-square
method: If we let $Y↓n$ be $2↑e X↓n$, so that $Y↓n$ is a double-precision
number obtained by placing $e$ zeros to the right of the binary
representation of $X↓n$, then $Y↓{n+1}$ consists of precisely
the middle $2e$ digits of $Y↑{2}↓{n} + 2↑e Y↓n$! In other words, Coveyou's
method is almost identical to a somewhat degenerate double-precision
middle-square method, yet it is guaranteed to have a long period;
further evidence of its randomness is proved in exercise
3.3.4--25.

\yskip
Other generalizations of Eq.\ (1) also suggest themselves; for
example, we might try to extend the period length of the sequence.
The period of a linear congruential sequence is extremely long;
when $m$ is approximately the word size of the computer, we
usually get periods on the order of $10↑9$ or more, so that typical
calculations will use only a very small portion of the sequence.
On the other hand, when we discuss the idea of ``accuracy'' in Section
3.3.4 we will see that the period length influences the degree of
randomness achievable in a sequence. Therefore it is occasionally desirable to
seek a longer period, and several methods are available for
this purpose. One technique is to make $X↓{n+1}$ depend on both $X↓n$
and $X↓{n-1}$, instead of just on $X↓n$; then the period
length can be as high as $m↑2$, since the sequence will not
begin to repeat until we have $(X↓{n+λ}, X↓{n+λ+1}) = (X↓n,
X↓{n+1}).$

The simplest sequence in which $X↓{n+1}$ depends
on more than one of the preceding values is the Fibonacci sequence,
$$X↓{n+1} = (X↓n + X↓{n-1})\mod m.\eqno (5)$$
This generator was considered in the early 1950s, and
it usually gives a period length greater than $m$; but tests
have shown that the numbers produced by the Fibonacci recurrence
(5) are definitely {\sl not} satisfactorily random, and so at
the present time the main interest in (5) as a source of random
numbers is that it makes a nice ``bad example.'' We may also
consider generators of the form
$$X↓{n+1} = (X↓n + X↓{n-k})\mod m,\eqno (6)$$
when $k$ is a comparatively large value. These were introduced by
Green, Smith, and Klem [{\sl JACM \bf 6} (1959), 527--537],
who reported that, when $k ≤ 15$, the sequence fails to pass
the ``gap test'' described in Section 3.3.2, although when $k
= 16$ the test was satisfactory.

A much better type of additive generator was devised in 1958
by G. J. Mitchell and D. P. Moore [unpublished], who
suggested the somewhat unusual sequence defined by
$$X↓n = (X↓{n-24} + X↓{n-55})\mod m,\qquad n ≥ 55,\eqno(7)$$
where $m$ is even, and where $X↓0, \ldotss,
X↓{54}$ are arbitrary integers not all even. The constants
24 and 55 in this definition were not chosen at random, they
are special values that happen to have the property that the
least significant bits $\langle X↓n \mod 2 \rangle$ will have
a period of length $2↑{55} - 1$. Therefore the sequence $\langle
X↓n\rangle$ must have a period at least this long. Exercise
11, which explains how to calculate the period length of such
sequences, proves that (7) has a period of length $2↑f(2↑{55}
- 1)$ for some $f,\ 0 ≤ f < e$, when $m = 2↑e$.

At first glance Eq.\ (7) may not seem to be extremely well suited
to machine implementation, but in fact there is a very efficient
way to generate the sequence using a cyclic list:

\algbegin Algorithm A (Additive number generator). Initially
locations $Y[1], Y[2], \ldotss,
Y[55]$ are set to the values $X↓{54}, X↓{53}, \ldotss, X↓0$, respectively;
$j$ is set to 24 and $k$ is set to 55. Successive performances
of this algorithm will produce the numbers $X↓{55}, X↓{56}, \ldots$
as output.

\algstep A1. [Add.] (If we are about to
output $X↓n$ at this point, $Y[j]$ now equals $X↓{n-24}$
and $Y[k]$ equals $X↓{n-55}$.)\xskip Set $Y[k] ← (Y[k] +
Y[j])\mod 2↑e$, and output $Y[k]$.

\algstep A2. [Advance.] Decrease $j$ and $k$ by 1. If now
$j = 0$, set $j ← 55$; otherwise if $k = 0$, set $k ← 55.$\quad\blackslug

\yyskip\noindent This algorithm in \MIX\ is simply the following:

\algbegin Program A (Additive number generator). Assuming that index
registers 5 and 6 are not touched by the remainder of the program in
which this routine is embedded, the following code performs
Algorithm A and leaves the result in register A\null. $\rI5 ≡ j$, $\rI6 ≡ k.$
$$\vcenter{\mixtwo{
LDA⊗Y,6⊗\understep{A1. Add.}\cr
ADD⊗Y,5⊗$Y↓k+Y↓j$ (overflow possible)\cr
STA⊗Y,6⊗\qquad$→Y↓k$.\cr
DEC5⊗1⊗\understep{A2. Advance.} $j ← j - 1$.\cr
DEC6⊗1⊗$k ← k - 1$.\cr
J5P⊗*+2\cr
ENT5⊗55⊗If $j = 0$, set $j ← 55$.\cr
J6P⊗*+2\cr
ENT6⊗55⊗If $k = 0$, set $k ← 55$.\quad\blackslug\cr}}$$

\yyskip This generator is usually faster than
the previous methods we have been discussing, since it does
not require any multiplication. Besides its speed, it has the
longest period we have seen yet; and it has consistently produced
reliable results, in extensive tests since its invention in
1958. Furthermore, as Richard Brent has observed, it can be made to
work correctly with floating-point numbers, avoiding the need to
convert between integers and fractions (cf.\ exercise 23).
Therefore it may well prove to be the very {\sl best}
source of random numbers for practical purposes. The only reason it is
difficult to recommend sequence (7) wholeheartedly is that there is
still very little theory to prove that it does or does not have
desirable randomness properties; essentially all we know for sure is that
the period is very long, and this is not enough. John Reiser (Ph.\ D. thesis,
Stanford Univ., 1977) has shown, however, that an additive sequence like
(7) will be well distributed in high dimensions, provided that a certain
plausible conjecture is true (cf.\ exercise 26).
%folio 37 galley 11 (C) Addison-Wesley 1978	*
The fact that the special numbers (24, 55)
in (7) work so well follows from theory developed in some of
the exercises below. Table 1 lists all pairs $(l, k)$ for which
the sequence $X↓n = (X↓{n-l} + X↓{n-k})\mod 2$ has period length
$2↑k - 1$, when $k < 100$. The pairs $(l, k)$ for small as well as larger
$k$ are shown, for the sake of completeness; the
pair (1, 2) corresponds to the Fibonacci sequence mod 2, whose
period has length 3. However, only pairs $(l, k)$ for relatively large
$k$ should be used to generate random numbers in practice.

\topinsert{\tablehead{Table 1}
\ctrline{\hbox expand 9pt{SUBSCRIPT PAIRS YIELDING LONG PERIODS MOD 2}}
\vskip 3 pt plus 2 pt minus 2 pt
\hrule
\vskip 3 pt plus .5 pt minus 1 pt
\tabskip 0pt plus 1pt
\halign to size{\lft{#}⊗\rt{#}⊗\rt{#}⊗\rt{#}⊗\rt{#}⊗\rt{#}⊗
\rt{#}⊗\rt{#}⊗\rt{#}\cr
(1, 2)⊗(1, 15)⊗(5, 23)⊗(7, 31)⊗(5, 47)⊗(21, 52)⊗(18, 65)⊗(28, 73)⊗(2, 93)\cr
(1, 3)⊗(4, 15)⊗(9, 23)⊗(13, 31)⊗(14, 47)⊗(24, 55)⊗(32, 65)⊗(31, 73)⊗(21, 94)\cr
(1, 4)⊗(7, 15)⊗(3, 25)⊗(13, 33)⊗(20, 47)⊗(7, 57)⊗(9, 68)⊗(9, 79)⊗(11, 95)\cr
(2, 5)⊗(3, 17)⊗(7, 25)⊗(2, 35)⊗(21, 47)⊗(22, 57)⊗(33, 68)⊗(19, 79)⊗(17, 95)\cr
(1, 6)⊗(5, 17)⊗(3, 28)⊗(11, 36)⊗(9, 49)⊗(19, 58)⊗(6, 71)⊗(4, 81)⊗(6, 97)\cr
(1, 7)⊗(6, 17)⊗(9, 28)⊗(4, 39)⊗(12, 49)⊗(1, 60)⊗(9, 71)⊗(16, 81)⊗(12, 97)\cr
(3, 7)⊗(7, 18)⊗(13, 28)⊗(8, 39)⊗(15, 49)⊗(11, 60)⊗(18, 71)⊗(35, 81)⊗(33, 97)\cr
(4, 9)⊗(3, 20)⊗(2, 29)⊗(14, 39)⊗(22, 49)⊗(1, 63)⊗(20, 71)⊗(13, 84)⊗(34, 97)\cr
(3, 10)⊗(2, 21)⊗(3, 31)⊗(3, 41)⊗(3, 52)⊗(5, 63)⊗(35, 71)⊗(13, 87)⊗(11, 98)\cr
(2, 11)⊗(1, 22)⊗(6, 31)⊗(20, 41)⊗(19, 52)⊗(31, 63)⊗(25, 73)⊗(38, 89)⊗(27, 98)\cr}
\vskip 3 pt plus .5 pt minus 1 pt
\hrule
\vskip 3 pt plus .5 pt minus 1 pt
\eightpoint
\hbox to size{For each pair $(l, k)$, the pair $(k - l, k)$ is also
valid (see exercise 24), hence only values of $l ≤ k/2$
are listed here. For extensions of this table, see N. Zierler
and J. Brillhart, {\sl Information and Control \bf 13} (1968),
541--554; {\bf 14} (1969), 566--569; {\bf 15} (1969), 67--69.}}

\yskip
Instead of considering only additive sequences, we can construct
useful random-number generators by taking general linear combinations
of $X↓{n-1}$, $\ldotss$, $X↓{n-k}$ for small $k$. In this case the
best results occur when the modulus $m$ is a large {\sl prime}\/; for example, $m$
can be chosen to be the largest prime number that
fits in a single computer word (see Table 4.5.4--1). When $m
= p$ is prime, the theory of finite fields tells us that it
is possible to find multipliers $a↓1, \ldotss, a↓k$ such that
the sequence defined by
$$X↓n = (a↓1X↓{n-1} + \cdots + a↓kX↓{n-k})\mod p\eqno(8)$$
has period length $p↑k - 1$; here $X↓0, \ldotss,
X↓{k-1}$ may be chosen arbitrarily but not all zero.\xskip (The
special case $k = 1$ corresponds to a multiplicative congruential
sequence with prime modulus, with which we are already familiar.)\xskip
The constants $a↓1, \ldotss, a↓k$ in (8) have the desired property
if and only if the polynomial
$$f(x) = x↑k - a↓1x↑{k-1} - \cdots - a↓k\eqno(9)$$
is a ``primitive polynomial modulo $p$,''
that is, if and only if this polynomial
has a root that is a primitive element of the field
with $p↑k$ elements (see exercise \hbox{4.6.2--16}).

Of course, the mere fact that suitable constants $a↓1$, $\ldotss$,
$a↓k$ {\sl exist} giving a period of length $p↑k - 1$
is not enough for practical purposes; we must be able to {\sl find}
them, and we can't simply try all $p↑k$ possibilities, since
$p$ is on the order of the computer's word size. Fortunately
there are exactly $\varphi (p↑k - 1)/k$ suitable choices of
$(a↓1, \ldotss, a↓k)$, so there is a fairly good chance of hitting
one after making a few random tries. But we also need a way
to tell quickly whether or not (9) is a primitive polynomial
modulo $p$; it is certainly unthinkable to generate up to $p↑k
- 1$ elements of the sequence and wait for a repetition! Methods
of testing for primitivity modulo $p$ are discussed by Alanen
and Knuth in {\sl Sankhy\=a} series A, {\bf 26} (1964),
305--328; the following criteria can be used: Let $r = (p↑k
- 1)/(p - 1)$.

\yskip\hang\textindent{i)}$(-1)↑{k-1}a↓k$ must be a primitive
root modulo $p$.\xskip (Cf.\ Section 3.2.1.2.)

\hang\textindent{ii)}The polynomial $x↑r$ must be congruent to $(-1)↑{k-1}a↓k$,
modulo $f(x)$ and $p$.

\hang\textindent{iii)}The degree of $x↑{r/q} \mod f(x)$, using polynomial arithmetic
modulo $p$, must be positive, for each prime divisor $q$ of $r$.

\yskip Efficient ways to compute the polynomial $x↑n \mod f(x)$,
using polynomial arithmetic modulo a given prime $p$, are discussed
in Section 4.6.2.

In order to carry out this test, we need to know the prime factorization
of $r = (p↑k - 1)/(p - 1)$, and this is the limiting factor
in the calculation; $r$ can be factored in a reasonable amount
of time when $k = 2$, 3, and perhaps 4, but higher values of
$k$ are difficult to handle when $p$ is large. Even $k = 2$
essentially doubles the number of ``significant random digits''
over what is achievable with $k = 1$, so larger values of $k$
will rarely be necessary.

An adaptation of the spectral test (Section 3.3.4) can be used
to rate the sequence of numbers generated by (8); see exercise
3.3.4--26. The considerations of that section show that we should
{\sl not} make the obvious choice of $a↓1 = +1$ or $-1$ when
it is possible to do so; it is better to pick large, essentially
``random,'' values of $a↓1, \ldotss, a↓k$ that satisfy the conditions,
and to verify the choice by applying the spectral test. A significant
amount of computation is involved in finding $a↓1, \ldotss, a↓k$,
but all known evidence indicates that the result will be a very
satisfactory source of random numbers. We essentially achieve
the randomness of a linear congruential generator with $k$-tuple
precision, using only single precision operations.

The special case $p = 2$ is of independent interest. Sometimes
a random-number generator is desired that merely produces a
random sequence of {\sl bits}---zeros and ones---instead of
fractions between zero and one. There is a simple way to generate
a highly random bit sequence on a binary computer, manipulating
$k$-bit words: Start with an arbitrary nonzero binary word
\.X. To get the next random bit of the sequence, do the following
operations, shown in \MIX's language (see exercise 16):
$$\vcenter{\mixtwo{
LDA⊗X⊗(Assume that overflow is now ``off.'')\cr
ADD⊗X⊗Shift left one bit.\cr
JNOV⊗*+2⊗Jump if high bit was originally zero.\cr
XOR⊗A⊗Otherwise adjust number with ``exclusive or.''\cr
STA⊗X⊗\qquad\blackslug\cr}}\eqno(10)$$
The fourth instruction here is the ``exclusive
or'' operation found on nearly all binary computers (cf.\ exercise
2.5--28 and Section 7.1);
it changes each bit position in which location \.A has a ``1''
bit. The value in location \.A is the binary constant $(a↓1 \ldotsm
a↓k)↓2$, where $x↑k - a↓1x↑{k-1} - \cdots - a↓k$
is a primitive polynomial modulo 2 as above. After the code
(10) has been executed, the next bit of the generated sequence
may be taken as the least significant bit of \.X (or, alternatively,
we could consistently use the most significant bit of \.X, if it were
more convenient to do so).
%folio 42 galley 12 (C) Addison-Wesley 1978	*
For example, consider Fig. 1, which
illustrates the sequence generated for $k = 4$ and $\.{CONTENTS(A)}=
(0011)↓2$. This is, of course, an unusually small value for
$k$. The right-hand column shows the sequence of bits of the
sequence, namely $1101011110001001\ldotsm$, repeating in a period
of length $2↑k - 1 = 15$.
This sequence is quite random, considering
that it was generated with only four bits of memory; to see
this, consider the adjacent sets of four bits occurring in the
period, namely 1101, 1010, 0101, 1011, 0111, 1111, 1110, 1100,
1000, 0001, 0010, 0100, 1001, 0011, 0110. In general, every
possible adjacent set of $k$ bits occurs exactly once in the
period, except the set of all zeros, since the period length
is $2↑k - 1$; thus, adjacent sets of $k$ bits are essentially
independent. We shall see in Section 3.5 that this is a very
strong criterion for randomness when $k$ is, say, 30 or more. Theoretical
results illustrating the randomness of this sequence are given
in an article by R. C. Tausworthe, {\sl Math.\ Comp.\ \bf 19}
(1965), 201--209.

\topinsert{
\hbox to size{
\hskip 19pt
\vbox{\eightpoint\vskip 12pt\halign{#\cr
1011\cr 0101\cr 1010\cr 0111\cr 1110\cr 1111\cr 1101\cr 1001\cr
0001\cr 0010\cr 0100\cr 1000\cr 0011\cr 0110\cr 1100\cr 1011\cr}}
\hfill
\vbox{
\hbox to 94mm{\caption Fig.\ 1. Successive contents of the computer word
\.X in the binary method, assuming that $k=4$ and $\.{CONTENTS(A)}
=(0011)↓2$.}
\vskip 15pt}
}}	% end \topinsert

Primitive polynomials of degree $≤100$ modulo 2 have been tabulated by
E. J. Watson, {\sl Math.\ Comp.\ \bf 16} (1962), 368--369. When
$k = 35$, we may take
$$\.{CONTENTS(A)} = (00000000000000000000000000000000101)↓2,$$
but the considerations of exercises 18 and 3.3.4--26 imply that
it would be preferable to find essentially ``random'' constants
that define primitive polynomials modulo 2.

{\sl Caution:} Several people have been trapped into believing
that this random-bit-generation technique can be used to generate
random whole-word fractions $(.X↓0X↓1 \ldotsm X↓{k-1})↓2$, $(.X↓kX↓{k+1}
\ldotsm X↓{2k-1})↓2$, $\ldotss$; but it is actually a poor source
of random fractions, even though the bits are individually quite
random.\xskip (See exercise 18.)

The concept of primitive polynomial modulo 2 is what underlies
the Mitchell--Moore generator (7); $x↑{55} + x↑{24} + 1$ is primitive,
and Table 1 is essentially a listing of all the primitive trinomials
modulo 2. A generator almost identical to that of Mitchell and
Moore was independently discovered in 1971 by T. G. Lewis and
W. H. Payne [cf.\ {\sl JACM \bf 20} (1973), 456--468], but using
``exclusive or'' instead of addition so that the period is exactly
$2↑{55} - 1$; each bit position in their generated numbers runs
through the same periodic sequence, but has its own starting
point.

We have now seen that sequences with $0 ≤ X↓n < m$ and period
$m↑k - 1$ can be found, when $X↓n$ is a suitable function of
$X↓{n-1}, \ldotss, X↓{n-k}$ and when $m$ is prime. The highest
conceivable period for {\sl any} sequence defined by a relation
of the form
$$X↓n = f(X↓{n-1}, \ldotss, X↓{n-k}),\qquad 0 ≤ X↓n < m,\eqno(11)$$
is easily seen to be $m↑k$. M. H. Martin [{\sl
Bull.\ Amer.\ Math.\ Soc.\ \bf 40} (1934), 859--864] was the first
person to show that functions achieving this maximum period
are possible for all $m$ and $k$; his method is easy to state,
but it is unfortunately not suitable for programming (see exercise
17). From a computational standpoint, the simplest known functions
$f$ that yield the maximum period $m↑k$ appear in exercise
21; the corresponding programs are, in general, not as efficient
for random-number generation as other methods we have described,
but they do give demonstrable randomness when the period as
a whole is considered.

\yskip Another important class of techniques deals with the {\sl combination}
of random-number generators, to get ``more random'' sequences.
There will always be people who feel that the linear congruential
methods, additive methods, etc., are all too simple to give
sufficiently random sequences; and it may never be possible to
{\sl prove} that their skepticism is unjustified (although we
believe it is), so it is pretty useless to argue the point.
There are reasonably efficient methods for combining two sequences
into a third one that should be haphazard enough to satisfy all but the most
hardened skeptic.

Suppose we have two sequences $X↓0, X↓1, \ldots$ and $Y↓0,
Y↓1, \ldots$ of random numbers between 0 and $m - 1$, preferably
generated by two unrelated methods. One suggestion has been
to add them together, mod $m$, obtaining the sequence $Z↓n =
(X↓n + Y↓n)\mod m$; in this case, the period will be quite long if the
period lengths of $\langle X↓n\rangle$ and $\langle Y↓n\rangle$ are
relatively prime to each other (see exercises 13 and 3.2.1.2--4).
Another approach, based on circular shifting and ``exclusive or''
on a binary computer, has been suggested by W. J. Westlake,
{\sl JACM \bf 14} (1967), 337--340.

A considerably different method has been suggested by M. D. MacLaren
and G. Marsaglia [{\sl JACM \bf 12} (1965), 83--89; {\sl CACM \bf 11}
(1968), 759], who use one random sequence to permute the elements
of another:

\algbegin Algorithm M (Randomizing by shuffling). 
Given methods for generating two sequences $\langle X↓n\rangle$
and $\langle Y↓n\rangle $, this algorithm will successively
output the terms of a ``considerably more random'' sequence.
We use an auxiliary table $V[0]$,
$V[1]$, $\ldotss$, $V[k - 1]$, where $k$ is some number chosen for
convenience, usually in the neighborhood of 100. Initially,
the $V$-table is filled with the first $k$ values of the $X$-sequence.

\algstep M1. [Generate $X, Y$.] Set $X$ and $Y$
equal to the next members of the sequences $\langle X↓n\rangle$
and $\langle Y↓n\rangle $, respectively.

\algstep M2. [Extract $j$.] Set $j ← \lfloor kY/m\rfloor $,
where $m$ is the modulus used in the sequence $\langle Y↓n\rangle$;
that is, $j$ is a random value, $0 ≤ j < k$, determined by $Y$.

\algstep M3. [Exchange.] Output $V[j]$ and then set $V[j] ← X$.\quad\blackslug

\yyskip 
As an example, assume that Algorithm
M is applied to the following two sequences, with $k = 64$:
$$\eqalign{X↓0 = 5772156649,⊗\qquad X↓{n+1} = (3141592653X↓n
+ 2718281829)\mod 2↑{35};\cr
Y↓0 = 1781072418,⊗\qquad Y↓{n+1} = (2718281829Y↓n + 3141592653)\mod 2↑{35}.\cr}
\eqno(12)$$
On intuitive grounds it appears safe to predict that
the sequence obtained by applying Algorithm M will satisfy virtually
{\sl anyone's} requirements for randomness in a computer-generated
sequence, because the relationship between nearby terms of the
output has been almost entirely obliterated. Furthermore, the
time required to generate this sequence is only slightly more
than twice as long as it takes to generate the sequence $\langle
X↓n\rangle$ alone.

Exercise 15 proves that the period length of Algorithm M's output
will be the least common multiple of the period lengths of $\langle
X↓n\rangle$ and $\langle Y↓n\rangle $, in most situations of
practical interest. In particular, the above example will have
a period of length $2↑{35}$.

However, there is an even better way to shuffle the elements of a sequence,
discovered by Carter Bays and S. D. Durham [{\sl ACM Trans.\ Math.\ 
Software \bf 2} (1976), 59-64]. Their approach, although it appears to be
superficially similar to Algorithm M\null, can give surprisingly better
performance even though it requires only one input sequence $\langle X↓n \rangle$
instead of two:

\algbegin Algorithm B (Randomizing by shuffling). Given a method for
generating a sequence $\langle X↓n \rangle$, this algorithm will successively
output the terms of a ``considerably more random'' sequence, using an
auxiliary table $V[0]$, $V[1]$, $\ldotss$, $V[k-1]$ as in Algorithm M\null.
Initially the $V$-table is filled with the first $k$ values of the $X$-sequence,
and an auxiliary variable $Y$ is set equal to the $(k+1)$st value.

\algstep B1. [Extract $j$.] Set $j←\lfloor kY/m\rfloor$, where $m$ is the modulus
used in the sequence $\langle X↓n\rangle$; that is, $j$ is a random value,
$0≤j<k$, determined by $Y$.

\algstep B2. [Exchange.] Set $Y←V[j]$, output $Y$, and then set $V[j]$ to the
next member of the sequence $\langle X↓n\rangle$.\quad\blackslug

\yyskip The reader is urged to work exercise 3, in order to get a feeling for
the difference between Algorithms M and B.

On \MIX\ we may implement Algorithm B by taking $k$ equal to the byte size,
obtaining the following simple generation scheme once the initialization
has been done:$$
\vcenter{\mixtwo{
LD6⊗Y(1:1)⊗$j←\hbox{high-order byte of }Y$.\cr
LDA⊗X⊗$\rA←X↓n$.\cr
INCA⊗1⊗(cf.\ exercise 3.2.1.1--1)\cr
MUL⊗A⊗$\rX←X↓{n+1}$.\cr
STX⊗X⊗``$n←n+1$.''\cr
LDA⊗V,6\cr
STA⊗Y⊗$Y←V[j]$.\cr
STX⊗V,6⊗$V[j]←X↓n$.\quad\blackslug\cr}}\eqno(13)$$

The output appears in register A\null. Note that Algorithm B requires
only four instructions of overhead per generated number.

F. Gebhardt [{\sl Math.\ Comp.\ \bf 21} (1967), 708--709] found
that satisfactory random sequences were produced by Algorithm
M even when it was applied to a sequence as nonrandom as the
Fibonacci sequence, with $X↓n = F↓{2n} \mod m$ and $Y↓n = F↓{2n+1}
\mod m$. However, it is also possible for Algorithm M to produce
a sequence {\sl less} random than the original sequences, if
$\langle X↓n\rangle$ and $\langle Y↓n\rangle$ are strongly related,
as shown in exercise 3. Such problems do not seem to arise with
Algorithm B\null. Since Algorithm B won't make a sequence any
less random, and since it probably enhances the randomness
substantially with very little extra cost, it can be recommended for
use in combination with any other random-number generator.

\exbegin{EXERCISES}
\trexno 1. [12] In practice,
we form random numbers using $X↓{n+1} = (aX↓n +
c)\mod m$, where the $X$'s are {\sl integers}, afterwards treating
them as the {\sl fractions} $U↓n = X↓n/m$. The recurrence relation
for $U↓n$ is actually
$$U↓{n+1} = (aU↓n + c/m)\mod 1.$$ Discuss the generation
of random sequences using this relation {\sl directly}, by making
use of floating-point arithmetic on the computer.

\trexno 2. [M20] A good source of random numbers will have $X↓{n-1}
< X↓{n+1} < X↓n$ about one-sixth of the time, since each of
the six possible relative orders of $X↓{n-1}$, $X↓n$, and $X↓{n+1}$
should be equally probable. However, show that the above ordering
{\sl never} occurs if the Fibonacci sequence (5) is used.

\exno 3. [23] (a) What sequence comes from Algorithm M if
$$X↓0 = 0,\quad X↓{n+1} = (5X↓n + 3)\mod 8,\qquad
Y↓0 = 0,\quad Y↓{n+1} = (5Y↓n + 1)\mod 8,$$
and $k = 4$?\xskip (Note that the potency is two, so $\langle X↓n\rangle$
and $\langle Y↓n\rangle$ aren't extremely random to start with.)\xskip
(b) What happens if Algorithm B is applied to this same sequence
$\langle X↓n\rangle$ with $k=4$?

\exno 4. [00] Why is the most significant byte used in the first line
of program (13), instead of some other byte?

\trexno 5. [20] Discuss using $X↓n = Y↓n$ in Algorithm M\null, in order
to improve the speed of generation. Is the result analogous to Algorithm B?

\exno 6. [10] In the binary method (10), the text states that
the low-order bit of \.X is random, if the code is performed repeatedly.
Why isn't the entire {\sl word\/} \.X random?

\exno 7. [20] Show that the full sequence of length $2↑e$ (i.e., a sequence in
which each of the $2↑e$ possible sets of $e$ adjacent bits occurs
just once in the period) may be obtained if program (10) is
changed to the following:
$$\vbox{\mixtwo{
LDA⊗X\cr
JANZ⊗*+2\cr
LDA⊗A\cr
ADD⊗X\cr
JNOV⊗*+3\cr
JAZ⊗*+2\cr
XOR⊗A\cr
STA⊗X\hskip .5 in \blackslug\cr}}$$
%folio 46 galley 13 (C) Addison-Wesley 1978	*
\exno 8. [M39] Prove
that the quadratic congruential sequence (3) has period length
$m$ if and only if the following conditions are satisfied:

\yskip{\sl
\halign{\rt{\qquad\rm#\xskip}⊗\lft{#}\cr
i)⊗$c$ is relatively prime to $m$;\cr
ii)⊗$d$ and $a - 1$ are both multiples of $p$, for all odd primes
$p$ dividing $m$;\cr
iii)⊗$d$ is even, and $d ≡ a-1 \modulo{4}$, if $m$ is a multiple of $4$;\cr
⊗$d ≡ a-1 \modulo{2}$, if $m$ is a multiple of $2$;\cr
iv)⊗either $d ≡ 0$ or both $a ≡ 1$ and $cd ≡ 6 \modulo{9}$,
if $m$ is a multiple of $9$.\cr}}	% end \sl
\yskip\noindent [{\sl Hint:} The sequence defined
by $X↓0 = 0,\ X↓{n+1} = dX↑{2}↓{n} + aX↓n + c$ has period length
$m$, modulo $m$, only if its period length is $r$ modulo any
divisor $r$ of $m$.]

\trexno 9. [M24] (R. R. Coveyou.)\xskip
Use the result of exercise 8 to prove that the modified middle-square
method (4) has a period of length $2↑{e-2}$.

\exno 10. [M29] Show that if $X↓0$ and $X↓1$ are not both even
and if $m = 2↑e$, the period of the Fibonacci sequence (5) is
$3 \cdot 2↑{e-1}$.

\exno 11. [M36] The purpose of this exercise is to analyze certain
properties of integer sequences satisfying the recurrence relation
$$X↓n = a↓1X↓{n-1} + \cdots + a↓kX↓{n-k},\qquad n ≥ k;$$
if we can calculate the period length of
this sequence modulo $m = p↑e$, when $p$ is prime, the period
length with respect to an arbitrary modulus $m$ is the least
common multiple of the period lengths for the prime power factors
of $m$.

\yskip
\def\\#1{\penalty0\;\biglp\hbox{modulo }f(z)\hbox{ and }#1\bigrp}
\hang\textindent{a)}If $f(z)$, $a(z)$, $b(z)$ are
polynomials with integer coefficients, let us write $a(z) ≡
b(z)\\m$ if $a(z) = b(z) + f(z)u(z)
+ mv(z)$ for some polynomials $u(z)$, $v(z)$ with integer coefficients.
Prove that when $f(0) = 1$ and $p↑e > 2$, ``If $z↑λ ≡ 1\\{p↑e}$,
$z↑λ \neqv 1\\{p↑{e+1}}$, then $z↑{pλ}≡1\\{p↑{e+1}}$, $z↑{pλ}\neqv 1\\{p↑{e+2}}$.''

\hang\textindent{b)}Let $f(z) = 1 - a↓1z - \cdots
- a↓kz↑k$, and let
$$G(z) = 1/f(z) = A↓0 + A↓1z + A↓2z↑2 + \cdotss.$$
Let $λ(m)$ denote the period
length of $\langle A↓n \mod m\rangle$. Prove that $λ(m)$ is
the smallest positive integer $λ$ such that $z↑λ ≡ 1\\m$.

\hang\textindent{c)}Given that $p$ is
prime, $p↑e > 2$, and $λ(p↑e) ≠ λ(p↑{e+1})$, prove
that $λ(p↑{e+r}) = p↑rλ(p↑e)$ for all $r ≥ 0$.\xskip $\biglp$Thus,
to find the period length of the sequence $\langle A↓n \mod
2↑e\rangle$, we can compute $λ(4)$, $λ(8)$, $λ(16)$, $\ldots$ until
we find the smallest $e ≥ 3$ such that $λ(2↑e) ≠
λ(4)$; then the period length is determined mod $2↑e$ for all
$e$. Exercise 4.6.3--26 explains how to calculate $X↓n$ for
large $n$ in $O(\log n)$ operations.$\bigrp$

\hang\textindent{d)}Show that any sequence of integers satisfying
the recurrence stated at the beginning of this exercise has
the generating function $g(z)/f(z)$, for some polynomial $g(z)$
with integer coefficients.

\hang\textindent{e)}Given that the polynomials $f(z)$ and $g(z)$ in part
(d) are relatively prime modulo $p$ (cf.\ Section 4.6.1), prove
that the sequence $\langle X↓n \mod p↑e\rangle$ has exactly
the same period length as the special sequence $\langle A↓n
\mod p↑e\rangle$ in (b).\xskip (No longer period could be obtained
by {\sl any} choice of $X↓0, \ldotss, X↓{k-1}$, since the general sequence
is a linear combination of ``shifts'' of the special sequence.)\xskip
[{\sl Hint:} By exercise 4.6.2--22 (Hensel's lemma), there
exist polynomials such that $a(z)f(z) + b(z)g(z) ≡ 1 \modulo{p↑e}$.]

\trexno 12. [M28] Find
integers $X↓0$, $X↓1$, $a$, $b$, and $c$ such that the sequence
$$X↓{n+1} = (aX↓n + bX↓{n-1} + c)\mod 2↑e,\qquad n ≥ 1,$$
has the longest period length of all sequences
of this type.\xskip [{\sl Hint:} It follows that $X↓{n+2} = \biglp (a + 1)X↓{n+1}
+ (b - a)X↓n - bX↓{n-1}\bigrp \mod 2↑e$; see exercise 11(c).]

\exno 13. [M20] Let $\langle X↓n\rangle$ and $\langle Y↓n\rangle$
be sequences of integers mod $m$ with periods of lengths $λ↓1$
and $λ↓2$, and form the sequence $Z↓n = (X↓n + Y↓n)\mod m$.
Show that if $λ↓1$ and $λ↓2$ are relatively prime, the sequence
$\langle Z↓n\rangle$ has a period of length $λ↓1λ↓2$.

\exno 14. [M24] Let $X↓n$, $Y↓n$, $Z↓n$, $λ↓1$, $λ↓2$ be as in the previous
exercise. Suppose that $λ↓1 = 2↑{e↓2}3↑{e↓3}5↑{e↓5}\ldots$
is the prime factorization of $λ↓1$, and similarly
suppose that $λ↓2 = 2↑{f↓2}3↑{f↓3}5↑{f↓5}\ldotss\,$.
Let $g↓p = \left(\max(e↓p, f↓p)\hbox{ if }e↓p ≠ f↓p,\hbox{ otherwise }0\right)$,
and let $λ↓0 = 2↑{g↓2}3↑{g↓3}5↑{g↓5}\ldotss\,$.
Show that the period length $λ↑\prime$ of the sequence $\langle Z↓n\rangle$ is a
multiple of $λ↓0$, but it is a divisor of $λ = \hbox{lcm}(λ↓1,λ↓2)$.
In particular, $λ↑\prime = λ$ if $(e↓p ≠ f↓p$ or $e↓p=f↓p=0)$ for each prime $p$.

\exno 15. [M27] Let the sequence $\langle X↓n\rangle$ in Algorithm
M have period length $λ↓1$, and assume that all elements of
its period are distinct. Let $q↓n = \min\leftset r\relv r > 0 $ and $\lfloor
kY↓{n-r}/m\rfloor = \lfloor kY↓n/m\rfloor\rightset$.
Assume that $q↓n < {1\over 2}λ↓1$
for all $n ≥ n↓0$, and that the sequence $\langle q↓n\rangle$
has period length $λ↓2$. Let $λ$ be the least common multiple of $λ↓1$
and $λ↓2$. Prove that the output sequence $\langle Z↓n\rangle$
produced by Algorithm M has a period of length $λ$.

\trexno 16. [M28] Let \.{CONTENTS(A)} in method (10) be $(a↓1a↓2
\ldotsm a↓k)↓2$ in binary notation. Show that the generated sequence
of low-order bits $X↓0$, $X↓1$, $\ldots$ satisfies the relation
$$X↓n = (a↓1X↓{n-1} + a↓2X↓{n-2} + \cdots + a↓kX↓{n-k})\mod 2.$$
[This may be regarded as another way to define
the sequence, although the connection between this relation
and the efficient code (10) is not apparent at first glance!]

\exno 17. [M33] (M. H. Martin, 1934.)\xskip Let $m$ and $k$ be positive integers,
and let $X↓1 = X↓2 = \cdots =
X↓k = 0$. For all $n > 0$, set $X↓{n+k}$ equal to the largest nonnegative
value $y < m$ such that the $k$-tuple $(X↓{n+1}, \ldotss, X↓{n+k-1},
y)$ has not already occurred in the sequence; in other words,
$(X↓{n+1}, \ldotss, X↓{n+k-1}, y)$ must differ from $(X↓{r+1},
\ldotss, X↓{r+k})$ for $0 ≤ r < n$. In this way, each possible
$k$-tuple will occur at most once in the sequence. Eventually
the process will terminate, when we reach a value of $n$ such
that $(X↓{n+1}, \ldotss, X↓{n+k-1}, y)$ has already occurred
in the sequence for all nonnegative $y < m$. For example, if $m =
3$ and $k = 3$ the sequence is 00022212202112102012001110100,
and the process terminates at this point.\xskip (a) Prove that when
the sequence terminates, we have $X↓{n+1} = \cdots
= X↓{n+k-1} = 0.$\xskip 
(b) Prove that {\sl every} $k$-tuple $(a↓1, a↓2, \ldotss, a↓k)$
of elements with $0 ≤ a↓j < m$ occurs in the sequence; hence
the sequence terminates when $n = m↑k$.\xskip [{\sl Hint:} Prove that
the $k$-tuple $(a↓1, \ldotss, a↓s, 0, \ldotss, 0)$ appears, when $a↓s ≠
0$, by induction on $s$.]\xskip Note that if we now define
$f(X↓n, \ldotss, X↓{n+k-1}) = X↓{n+k}$ for $1 ≤ n ≤ m↑k$, setting
$X↓{m↑k+k} = 0$, we obtain a function of maximum possible
period.

\exno 18. [M22] Let $\langle X↓n\rangle$ be the sequence of
bits generated by method (10),
with $k = 35$ and $\.{CONTENTS(A)}=(00000000000000000000000000000000101)↓2$.
Let $U↓n$ be the binary fraction $(.X↓{nk}X↓{nk+1}\ldotsm X↓{nk+k-1})↓2$;
show that this sequence $\langle U↓n\rangle$ fails the serial
test on pairs (Section 3.3.2B) when $d = 8$.

\exno 19. [M41] For each prime $p$ specified in the first column
of Table 1 in Section 4.5.4, find suitable constants $a↓1, a↓2$
as suggested in the text, such that the period length of (8),
when $k = 2$, is $p↑2 - 1$.\xskip (See Eq.\ 3.3.4--39 for an example.)

\exno 20. [M40] Calculate constants suitable for use as \.{CONTENTS(A)} in
method (10), having approximately the same number of zeros as
ones, for $2 ≤ k ≤ 64$.

\exno 21. [M35] (D. Rees.)\xskip The text explains how to find functions
$f$ such that the sequence (11) has period length $m↑k - 1$,
provided that $m$ is prime and $X↓0, \ldotss, X↓{k-1}$ are not
all zero. Show that such functions can be modified to obtain
sequences of type (11) with period length $m↑k$, for {\sl all} integers
$m$.\xskip [{\sl Hints:} Consider Lemma 3.2.1.2Q\null, the trick of exercise
7, and sequences such as $\langle pX↓{2n} + X↓{2n+1}\rangle$.]

\trexno 22. [M24] The text restricts discussion of the extended
linear sequences (8) to the case that $m$ is prime. Prove that
reasonably long periods can also be obtained when $m$ is ``square-free,''
i.e., the product of distinct primes.\xskip (Examination of Table
3.2.1.1--1 shows that $m = w \pm 1$ often satisfies this hypothesis;
many of the results of the text can therefore be carried over
to that case, which is somewhat more convenient for calculation.)

\trexno 23. [20] Discuss the sequence defined by $X↓n = (X↓{n-55}
- X↓{n-24})\mod m$ as an alternative to (7).

\exno 24. [M20] Let $0 < k < m$. Prove that the sequence defined
by $X↓n = (X↓{n-m+k} + X↓{n-m})\mod 2$ has period length $2↑m - 1$ whenever
the sequence defined by $Y↓n = (Y↓{n-k} + Y↓{n-m})\mod 2$ does.

\exno 25. [26] Discuss the alternative to Program A that changes
all 55 entries of the $Y$ table every 55th time a random number
is required.

\exno 26. [M48] (J. F. Reiser.)\xskip Let $p$ be prime and let $k$ be a positive
integer.  Given integers $a↓1,\ldotss,a↓k$ and $x↓1,\ldotss,x↓k$, let
$λ↓α$ be the period of the sequence $\langle X↓n\rangle$ generated by the
recurrence
$$X↓n=x↓n\mod p↑α,\quad 0≤n<k;\qquad X↓n=(a↓1X↓{n-1}+\cdots+a↓kX↓{n-k})\mod
p↑α,\quad n≥k;$$
and let $N↓α$ be the number of 0's that occur in the period (i.e., the number
of indices $j$ such that $\mu↓α≤j<\mu↓α+λ↓α$ and $X↓j=0$).
Prove or disprove the following conjecture: There exists a constant $c$
(depending possibly on $p$ and $k$ and $a↓1,\ldotss,a↓k$) such that
$N↓α≤cp↑{α(k-2)/(k-1)}$ for all $α$ and all $x↓1,\ldotss,x↓k$.

[{\sl Notes:} Reiser has proved that if the recurrence has maximum
period length mod $p$ (i.e., if $λ↓1=p↑k-1$), and if the conjecture holds,
then the $k$-dimensional discrepancy of $\langle X↓n\rangle$ will be
$O(α↑k p↑{-α/(k-1)})$ as $α→∞$; thus an additive generator like (7) would
be well distributed in 55 dimensions, when $m=2↑e$ and the entire period
is considered.\xskip (See Section 3.3.4 for the definition of discrepancy in
$k$ dimensions.)\xskip
The conjecture is a very weak condition, for if $\langle X↓n
\rangle$ takes on each value about equally often and if $λ↓α = p↑{α-1}(p↑k-1)$,
the quantity $N↓α \approx (p↑k-1)/p$ does not grow at all as $α$ increases.
Reiser has verified the conjecture for $k=3$. On the other hand he has shown
that it is possible to find unusually bad starting values $x↓1,\ldotss,x↓k$
(depending on $α$) so that $N↓{2α}≥p↑α$, provided that $λ↓α=p↑{α-1}(p↑k-1)$
and $k≥3$ and $α$ is sufficiently large.]

\exno 27. [M30] Suppose Algorithm B is being applied to a sequence $\langle X↓n
\rangle$ of period length $λ$, where $λ\grgr k$. Show that for fixed $k$ and all
sufficiently large $λ$, the output of the sequence will eventually be periodic
with the same period length $λ$, unless $\langle X↓n\rangle$ isn't very random to
start with.\xskip
[{\sl Hint:} Find a pattern of consecutive values of $\lfloor k X↓n
/m\rfloor$ that causes Algorithm B to ``synchronize'' its subsequent behavior.]

\exno 28. [40] (A. G. Waterman.)\xskip Experiment with linear congruential
sequences with $m$ the square or cube of the computer word size,
while $a$ and $c$ are single-precision numbers.

\vfill\eject